      subroutine filterfft(nvtx, ijkvtx, ibar, jbar, kbar, ncyc, ibp2, jbp2, &
         kbp2, bzv, iwid, jwid, kwid, idx, idy, idz) 
!-----------------------------------------------
!   M o d u l e s 
!-----------------------------------------------
      USE vast_kind_param, ONLY:  double 
      use modify_com_M 
 
!...Translated by Pacific-Sierra Research 77to90  4.3E  14:13:36   8/20/02  
!...Switches: -yf -x1             
      implicit none
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      integer  :: nvtx 
      integer , intent(in) :: ibar 
      integer , intent(in) :: jbar 
      integer  :: kbar 
      integer  :: ncyc 
      integer , intent(in) :: ibp2 
      integer , intent(in) :: jbp2 
      integer , intent(in) :: kbp2 
      integer , intent(in) :: iwid 
      integer , intent(in) :: jwid 
      integer , intent(in) :: kwid 
      integer  :: idx 
      integer , intent(in) :: idy 
      integer  :: idz 
      integer  :: ijkvtx(*) 
      real(double) , intent(inout) :: bzv(*) 
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      integer :: ibp1, jbp1, kbp1, i, j, k, ijk, ngx, ngy, kx, ky 
      real(double) :: pi, dx, dy, factor, error, dnorm 
      complex, dimension(idx,idy) :: ffta, fftc, fftb 
!-----------------------------------------------
 
      ibp1 = ibp2 - 1 
      jbp1 = jbp2 - 1 
      kbp1 = kbp2 - 1 
 
      pi = 3.1415926536 
      dx = 2.*pi 
      dy = dx 
 
      do i = 2, ibp2 
         do j = 2, jbp2 
            ffta(i-1,j-1) = 0. 
            ffta(i-1,j-1) = ffta(i-1,j-1) + sum(bzv(kwid+1+(i-1)*iwid+(j-1)*&
               jwid:(kbp2-1)*kwid+1+(i-1)*iwid+(j-1)*jwid:kwid)) 
            ffta(i-1,j-1) = ffta(i-1,j-1)/float(kbp1) 
         end do 
      end do 
 
      call fft2d (ibp1, jbp1, ffta, idx, fftc, idx) 
 
      ngx = ibp1/2. 
      ngy = jbp1/2. 
 
      do i = 1, ibp1 
 
         if (i <= ngx) then 
            kx = (i - 1)/dx*2.*pi 
         else 
            kx = -(ibp1 - i + 1)/dx*2.*pi 
         endif 
 
         do j = 1, jbp1 
 
            if (j <= ngy) then 
               ky = (j - 1)/dy*2.*pi 
            else 
               ky = -(jbp1 - j + 1)/dy*2.*pi 
            endif 
            if (abs(kx)<=ibar/2 .and. abs(ky)<=jbar/2) cycle  
 
            fftc(i,j) = 0. 
 
         end do 
 
      end do 
 
      fftc(1,1) = 0. 
 
      call fft2b (ibp1, jbp1, fftc, idx, fftb, idx) 
 
      factor = 1./(ibp1*jbp1) 
 
      error = 0. 
      dnorm = 0. 
      do i = 1, ibp1 
         do j = 1, jbp1 
            error = error + abs(ffta(i,j)-fftb(i,j)*factor) 
            dnorm = dnorm + abs(ffta(i,j)) 
         end do 
      end do 
!      write(*,*)'test filter:',error/dnorm
 
      do i = 2, ibp2 
         do j = 2, jbp2 
            do k = 2, kbp2 
               ijk = 1 + (i - 1)*iwid + (j - 1)*jwid + (k - 1)*kwid 
               bzv(ijk) = dble(fftb(i-1,j-1))*factor 
            end do 
         end do 
      end do 
      return  
      end subroutine filterfft 
